clip-02-05

clip-02-05.jl

Load Julia packages (libraries) needed for the snippets in chapter 0

using StatisticalRethinking, Optim
gr(size=(600,300));
Plots.GRBackend()

snippet 3.2

p_grid = range(0, step=0.001, stop=1)
prior = ones(length(p_grid))
likelihood = [pdf(Binomial(9, p), 6) for p in p_grid]
posterior = likelihood .* prior
posterior = posterior / sum(posterior)
samples = sample(p_grid, Weights(posterior), length(p_grid));
samples[1:5]
5-element Array{Float64,1}:
 0.573
 0.836
 0.788
 0.634
 0.727

snippet 3.3

Draw 10000 samples from this posterior distribution

N = 10000
samples = sample(p_grid, Weights(posterior), N);
10000-element Array{Float64,1}:
 0.855
 0.604
 0.769
 0.653
 0.733
 0.683
 0.851
 0.721
 0.692
 0.751
 ⋮
 0.661
 0.548
 0.673
 0.636
 0.581
 0.753
 0.76
 0.46
 0.642

In StatisticalRethinkingJulia samples will always be stored in an MCMCChains.Chains object.

chn = MCMCChains.Chains(reshape(samples, N, 1, 1), ["toss"]);
Object of type Chains, with data of type 10000×1×1 Array{Float64,3}

Log evidence      = 0.0
Iterations        = 1:10000
Thinning interval = 1
Chains            = 1
Samples per chain = 10000
parameters        = toss

parameters
      Mean    SD   Naive SE  MCSE   ESS
toss 0.6352 0.1405   0.0014 0.0014 10000

Describe the chain

describe(chn)
Log evidence      = 0.0
Iterations        = 1:10000
Thinning interval = 1
Chains            = 1
Samples per chain = 10000
parameters        = toss

Empirical Posterior Estimates
────────────────────────────────────────
parameters
      Mean    SD   Naive SE  MCSE   ESS
toss 0.6352 0.1405   0.0014 0.0014 10000

Quantiles
────────────────────────────────────────
parameters
      2.5% 25.0% 50.0% 75.0% 97.5%
toss 0.098  0.54 0.644 0.741 0.977

Plot the chain

plot(chn)
0 2500 5000 7500 10000 0.2 0.4 0.6 0.8 1.0 toss Iteration Sample value 0.00 0.25 0.50 0.75 1.00 0.0 0.5 1.0 1.5 2.0 2.5 toss Sample value Density

snippet 3.4

Create a vector to hold the plots so we can later combine them

p = Vector{Plots.Plot{Plots.GRBackend}}(undef, 2)
p[1] = scatter(1:N, samples, markersize = 2, ylim=(0.0, 1.3), lab="Draws")
0 2500 5000 7500 10000 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Draws

snippet 3.5

Analytical calculation

w = 6
n = 9
x = 0:0.01:1
p[2] = density(samples, ylim=(0.0, 5.0), lab="Sample density")
p[2] = plot!( x, pdf.(Beta( w+1 , n-w+1 ) , x ), lab="Conjugate solution")
0.00 0.25 0.50 0.75 1.00 0 1 2 3 4 5 Sample density Conjugate solution

Add quadratic approximation

plot(p..., layout=(1, 2))
0 2500 5000 7500 10000 0.0 0.2 0.4 0.6 0.8 1.0 1.2 Draws 0.00 0.25 0.50 0.75 1.00 0 1 2 3 4 5 Sample density Conjugate solution

End of 03/clip-02-05.jl

This page was generated using Literate.jl.